Optimized Verlet-like algorithms for molecular dynamics simulations 
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New explicit velocity- and position- Verlet-like algorithms of the second order are proposed to 
integrate the equations of motion in many-body systems. The algorithms are derived on the basis 
of an extended decomposition scheme at the presence of a free parameter. The nonzero value for 
this parameter is obtained by reducing the influence of truncated terms to a minimum. As a result, 
the new algorithms appear to be more efficient than the original Verlet versions which correspond 
to a particular case when the introduced parameter is equal to zero. Like the original versions, the 
proposed counterparts are symplectic and time reversible, but lead to an improved accuracy in the 
generated solutions at the same overall computational costs. The advantages of the new algorithms 
are demonstrated in molecular dynamics simulations of a Lennard- Jones fluid. 
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The method of molecular dynamics (MD) is a powerful 
tool for the prediction and study of various phenomena 
in physics, chemistry and biology. In MD simulations we 
deal with the necessity to solve numerically the equations 
of motion for a many-body system composed of interact- 
ing particles. The most of traditional algorithms, such as 
Runge-Kutta and predictor-corrector schemes are 
usually unsuitable for integration of the resulting differ- 
ential equations, because the solutions obtained exhibit 
a high instability on MD scales of time ^ . 

A variety of alternative algorithms were proposed and 
implemented over the years |-§. These include the well- 
known velocity- Verlet (W) integrator This second- 
order integrator is employed in the great majority of MD 
simulations due to its simplicity and exceptional stabil- 
ity. Moreover, the VV algorithm is symplectic, time re- 
versible, and able to reach a high level of accuracy with 
minimal number of force evaluations per time step . 
In addition, the VV approach can be modified to inte- 
grate not only translational motion in atomic systems, 
but also simulate more complicated molecular and spin 
liquids §-11 . 

The question of how to improve the efficiency of inte- 
gration for atomic systems with long-range interactions 
has also been considered. As a result, so-called multi- 
ple time scale integrators have been introduced p3| , p^ . 
In these integrators, the additional slow subdynamics is 
treated in a specific way using the weakness of the long- 
range forces. The faster motion, caused by the interac- 
tions at short interparticle distances, remains to be inte- 
grated with the help of usual basic algorithms, such as 
VV integrator, for instance. 

In the present paper we show that even within the 
basic consideration of translational motion (when ad- 
ditional splitting of interaction potentials into multiple 
scale components is not longer allowed), the VV algo- 
rithm presents, in fact, only a particular case among a 
whole family of symplectic reversible integrators of the 



second-order. This case, appears to be not so optimal 
and more efficient second-order algorithms are possible. 

The equations of motion for a classical system consist- 
ing of N particles can be cast in the following compact 
form 



(1) 



where p= {r^, v^} denotes the full set (i = 1, 2, . . . , TV) of 
phase variables with and being the position and ve- 
locity, respectively, of the ith particle, L is the Liouville 
operator. 
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fi — ^J2j{j=£i) f' /'''ij designates the force acting 
on the particles of mass m each, due to the interactions 
described by the potential ip(rij), and = — Vj. The 
Liouville operator has been split in Eq. (2) into the free- 
motion A = v-d/dr and potential B = i/rrt-d/dv parts 
with V = {vj, r iE {r J, and f = {fj. 
The formal solution of Eq. (1) is 
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where h denotes the time step. Of course, the exponen- 
tial propagator e^^ cannot be evaluated exactly at any 
h (solutions in quadratures are possible only for = 2 
that is not relevant for our consideration of many-body 
systems when A^ 1). However, at small enough val- 
ues of h, the total propagator allows to be decomposed 
iT5|-ra as 
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where the coefScients ap and bp are chosen in such a way 
to provide the highest possible value for > 1 at a given 
integer number P > 1. Then, starting from an initial con- 
figuration p(0), the evolution of the system can be inves- 
tigated during arbitrary times t by repeating the single- 



step propagation, p{t) 
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(5) 



where I — t/h is the total number of steps and the trun- 
cation terms 0{h^^^), appearing in Eq. (4), have been 
neglected. 

The main advantage of decomposition (4) is that the 
exponential subpropagators e^"^ and e^"^ are analytically 
integrable. Indeed, 



e^V = e^"^/^''^{r, v} = {r + rv, v} , 
e^^p = ef/™-9/9vr|j.^ ^1 = {r, V + rf /m} 



(6) 



that represents simple shifts of positions and velocities, 
respectively (with t being equal to Uph or bph). In addi- 
tion, the generated trajectory (5) behaves symplectically 
(like exact solutions), because such separate shifts do not 
change the volume in phase space. The time reversibility 
S{~t)p{t) ~ p{0) of solutions (following from the prop- 
erty S~^{t) = S{—t) of evolution operators) can also be 
obtained by imposing additional conditions on the co- 
efficients ftp and bp. Namely, ai = 0, flp+i = ap-p+i. 



and bp = 6p_p+i, or Op 



ap-p-^i, and bp = ap^p at 



6p = 0. In other words, the subpropagators e^'^ and 
e^'^ should enter symmetrically in the decompositions. 
Then the even-order truncation terms 0{h?^) will dis- 
appear automatically in Eq. (4) at fc < K/2. For this 
reason, the order K of reversible algorithms may accept 
only even numbers. The cancellation of odd-order terms 
0{h?^^^) will be provided by satisfying a set of basic 
conditions for ap and bp. For example, the condition 

S^=i = S^=i &p = 1 is required to cancel the first- 
order truncation uncertainties. 

The method just highlighted is quite general to build 
numerical integrators of arbitrary orders. In particular, 
the second-order {K — 2) VV algorithm 



(7) 



is immediately reproduced from Eq. (4) at P = 2 and 
fli = 0, fei = ^2 = 1/2: 0'2 = 1- The case when the opera- 
tors A and B are replaced by each other [A ^ B) is also 
possible, and we come to the so-called position- Verlet 
(PV) algorithm [|l|, e('4+^)'' = e^'^/^e^'^e'^''/^ + ^^(/iS), 
corresponding to the choice oi = a2 = 1/2, bi = 1, 
and &2 = 0. Algorithms of higher orders can also be 
derived in a similar way. For instance, the fourth-order 
[K = 4) algorithm by Forest and Ruth |l^ is obtained 
from Eq. (4) at P = 4, whereas sixth-order {K = 6) 



schemes are derivable beginning from P = 8. The 
high-order schemes involve, however, too large number 
of force recalculations, and appear to be less efficient in 
most of MD applications than second-order algorithms. 

Despite the fact that the method of construction of 
time-reversible integrators using syniplectic decomposi- 
tions is not new, some important cases have never been 
considered and have been completely ignored in the lit- 
erature. This concerns, in particular, the following ques- 
tion. Are the above Verlet algorithms optimal in view 
of the time efficiency among all possible basic (i.e. with 
single splitting of the Liouville operator) decomposition 
integrators of the second-order? We can say only that 
the Verlet algorithms do minimize the number of force 
evaluations per time step. However, as will be shown 
below, this does not guarantee the optimization with re- 
spect to the overall number of force recalculations (the 
most time-consuming part of MD simulations), which are 
necessary to perform during a fixed observation time in 
order to achieve a given precision in solutions. 

It can be seen readily that the Verlet algorithms (P — 
2) require only one (P— 1) force evaluations per time step 
/i, whereas the fourth- and higher-order schemes P > 4 
need in three or more such evaluations. Let us consider 
now the intermediate case P — 3 which leads to an ex- 
tended time-reversible propagation in the form 

^iA+B)h ^ ^Aih^B^^Ail^2i)h^B^^Aih _^ (jj^3 ^ (r^^j^i^ 

(8) 

following from Eq. (4) at ai = 03 = ^, 02 = 1 — 2^, 
bi = b2 = 1/2, and 63 = 0. Again, the propagation with 
yl P is also acceptable (then ai = 0, 61 = 63 = ^, 
62 = 1 — 2f and 02 = 03 = 1/2). Formula (8) represents 
a whole family of symplectic time-reversible integrators 
of the second-order in which a particular member can 
be extracted by choosing a value for free parameter ^. 
For ^ = 0, Eq. (8) reduces to the VV (see Eq. (7)) or 
PV (at A^ B) algorithm. The extended (when ^ ^ 0) 
propagation requires already two, instead of one, force re- 
calculation per time step. For this reason, we can come 
to an incorrect conclusion that such a propagation has 
no advantage over the Verlet algorithms. 

In order to prove that the above conclusion is indeed 
incorrect, let us analyze in more detail the influence of 
truncation errors Ch^ on the result. Expanding both the 
sides of Eq. (8) into Taylor's series with respect to h, one 
finds 



C = a{0[A, [B,A]]+m[B,[B.A]], 



where 
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and [ , ] denotes the commutator of two operators. The 
norm of C with respect to the third-order commutators 
[A, [B,A]] and [B,[B,A]] is 
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7(0 = V«'(0+/3'(0- 



(11) 



Then the norm of local uncertainties Cph^ appearing in 
phase trajectory p during a single-step propagation given 
by Eqs. (3) and (8) can be expressed in terms of 7 and 
h as g — ^h^ . During a whole integration over a fixed 
time interval t, the total number / of such single steps 
is proportional to (see Eq. (5)). As a result, the 
local third-order uncertainties will accumulate step by 
step leading at < /i to the second-order global errors 
r = gh,-^, i.e., 



(12) 



Extended propagation (8) can now be optimized with 
respect to ^ by minimizing the function 7(5). As can be 
verified readily, the minimum of 7(^) is achieved at ^ = 
where 



1 (2^326 + 36)^^^ 
*' ~ 2 12 



0.1931833275037836 
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(13) 



and consists 7(C) w 0.00855. On the other hand, the 
value 7(0) of 7 corresponding to the Verlet algorithms 
(when ^ = 0) is equal to 7(0) « 0.0932, i.e., it increases in 
7(0)/7(C) « 11 times. Remembering that the extended 
propagation requires two force evaluation per time step h, 
it should be performed with double step size 2h with re- 
spect to that of the Verlet algorithms, in order to provide 
the same number of total force recalculations during the 
integration over the same time interval. Therefore, the 
extended propagation will be more efficient if the follow- 
ing inequality r(^, 2h) < r(^ = 0, h) takes place. Taking 
into account Eq. (12), such an inequality can be rewrit- 
ten as 7(0)/7(^) > 4, and thus it is fulfilled completely 
in the optimization regime. In particular. 



r(C,2/i) 



0.367, 



(14) 



indicating that the optimized propagation, being applied 
even with double sizes of the time step, will reduce the 
global errors approximately in three times. 

In view of Eqs. (3), (6), and (8), more explicit ex- 
pressions for the single-step propagation of position and 
velocity from time t to t+h within the optimized W-like 
algorithm are: 
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whereas the optimized PV-like algorithm (when A 
in Eq. (8)) reads: 
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Vl = v{t) + ^{[rmh 

ri = r(t) + viV2 

V2 =vi + if [ri] (1-2^/1 
r{t + h) =ri +V2/1/2 
v(t + /i)=V2 + ^f[r(t + /i)]e/i, 



(16) 



where the parameter ^ should take its optimal value ( 
(see Eq. (13)). The algorithms are simple, require only 
slight modification with respect to the original Verlet ver- 
sions, and can easily be implemented in program codes. 

It is worth pointing out that the order of local errors 
0{h^) = Cph^ = C{r, v}ft"^ remains three in both posi- 
tion r{t + h) and velocity r(t + h) for both the optimized 
algorithms (15) and (16) (because the functions a(^) and 
f}{^) cannot tend to zero simultaneously at any ^). Note 
also that the minimization of third-order uncertainties 
Ch^ in Eq. (8) automatically minimizes the fourth-order 
truncation terms 0{h'^) which are connected with C by 
the relation 0(/i4) = ^ {CiA + B) + {A + B)C)h^ + 0{h^). 
Further optimization is also possible in specific cases. For 
instance, some MD applications are aimed exclusively 
at the investigation of structural properties of the sys- 
tem. Then the accuracy of determining particle posi- 
tions will play more important role than that of veloci- 
ties. In such a situation, it is quite natural to increase 
the precision in evaluation of r{t + h) by reducing the 
position part Crh^ = (a(^)Ci + P{0C2)rh^ of third- 
order truncation errors to zero, where Ci = [A, [B, A]] 
and C2 = [B, [i?,^]] (see Eq. (9)). This reduction can 
be realized, because (as can be shown using the explicit 
expressions for A and B) the operator C2 vanishes when 
acting on position, i.e. C2r = 0, whereas Cir ^ (as 
well as Civ ^ and C2V ^ 0). The influence of C\r 
can be reduced to zero also by choosing such S, at which 
a(^) =0. Among the two roots (1=F l/V3)/2 of equation 
a(^) = 0, the preference should be given to the first of 
them, (1 — l/\/3)/2, because it leads to a smaller value 
for PiO- Then replacing C by (1 - l/V3)/2 in Eq. (15), 
we come to a positionally optimized VV-like algorithm in 
which the positions will be generated up to the fourth- 
order truncation imcertainties 0{h*). 

Another useful application of the positionally opti- 
mized algorithm is the case of weakly interacting sys- 
tems, where the Liouville operator can be presented in 
the form L = A + eB with e <C 1. Then the operator 
[B, [i?,^]] = C2, which forms the third-order errors in 
velocity, will be proportional to and, thus, can be ne- 
glected. For the same reason, the corresponding fourth- 
order uncertainties i(C2(A-FS)-F(^4-S)C2))/i** will also 
behave like e^. In such a case, the positionally optimized 
algorithm can be considered as a quasi fourth-order inte- 
grator which, contrary to the usual fourth-order schemes, 
will require only two (instead of three) force evaluations 
per time step. 

In order to obtain a positionally optimized algorithm 
within the PV-like integration (16), it is necessary simply 
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to replace ^ by the root 1/6 of equation (3{^) = 0. Note 
that the values 1/6 « 0.167 and (1 - l/V^)/2 « 0.211 
are close enough to the optimal solution (13) which mini- 
mizes the total position- velocity uncertainties. Neverthe- 
less, the positionally optimized algorithms are not recom- 
mended to be used in general case when both the position 
and velocity must be evaluated with a maximal accuracy. 
In other words, in such partially optimized algorithms the 
increased precision in position evaluation is achieved at 
the expense of decreasing accuracy in determining the 
velocities. Indeed, 7(0)/7(l/6) = 7(0)/|a(l/6)| w 7 and 
7(0)/7((l - l/V3)/2) ^ 7(0)/|/3((l - 1/V3)/2)| « 8 that 
is less than the factor 7(0) /7(C) w 11 corresponding to 
optimal value (13). 

Our theoretical predictions were verified by testing 
the VV- and PV-like algorithms in MD simulations of a 
Lennard- Jones fluid (LJ). We considered a system com- 
posed of = 256 particles interacting through the LJ po- 
tential <I>(r) — 4u[((T/r)-^^ — (cr/r)^] in a basic cubic box 
of volume V = using periodic boundary conditions. 
The LJ potential was truncated at rc — L/2 fa 3.36ct and 
shifted to be zero at the truncation point to avoid the 
force singularities, i.e., (p(r) = <I>(r) — $(rc) at r < Tc 
and f(r) = otherwise. The simulations were carried 
out in a microcanonical ensemble at a reduced density 
of 71* = ^(T'^ = 0.845 and a reduced temperature of 
T* = k^T /u = 1.7. The equations of motion were in- 
tegrated with the help of Eqs. (15) and (16) in which 
the parameter ^, being constant within each run, varied 
from one run to another. All the runs started from an 
identical well equilibrated initial configuration p(0), and 
further continued I = 10 000 time steps. The precision 
of generated solutions was measured in terms of the rel- 
ative total energy fluctuations £ = {{E — {E)Y) /\{E)\, 

where E = \ Yl!i=i '^^'i^ + \ J2i!^j viUj) and ( ) denotes 
the microcanonical averaging. Note that in microcanon- 
ical ensembles the total energy is an integral of motion, 
E{t) = E{0), and the above fluctuations should be equal 
to zero if the equations of motion are solved exactly. So 
that in approximate MD integration, smaller values of 
£ will correspond to a more precise evaluation of phase 
trajectory. 

The total energy fluctuations obtained in the simula- 
tions at the end of the runs for four (flxed within each 
run) undimensional time steps, h* = h{u/ma'^y^'^ = 
0.01, 0.005, 0.0025, and 0.001, are shown in Fig. 1 as 
depending on free parameter ^. The subsets (a) and 
(b) of this figure correspond to the VV- and PV-like in- 
tegration, respectively. As can be seen, all the depen- 
dencies £(^, /i) have one minimum which locates at the 
same point ^ « 0.19 independently on the size h of the 
time step. This point coincides completely with the min- 
imum at C (Eq. (13)) of function 7(^) (Eq. (11)) which 
is included in Fig. 1 as well (dashed curves in the sub- 
sets). Moreover, the energy fluctuations £{^,h) appear 
to be proportional to the norm T{£^,h) of global errors 
(see Eq. (12)), and the coefficient of this proportionality 
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FIG. L The total energy fluctuations obtained in the simu- 
lations for different values of free parameter ^ at four reduced 
time steps, h* = 0.01, 0.005, 0.0025, and 0.001, using the VV- 
(subset (a)) and PV- (subset (b)) like integration (Eqs. (15) 
and (16), respectively). The simulation results are presented 
by circles connected by the solid curves. The function 7(^) 
(see Eq. (11)) is plotted in both the subsets by the dashed 
curve. 

almost does not depend on ^ and h. In addition, at each 
step size considered the energy fluctuations decrease at 
the minimum more than in ten times with respect to 
those at ^ = 0, that is in agreement with our predicted 
value 7(0)/7(C) w 11. 

The result for the total energy fluctuations as functions 
of the length of the simulations corresponding to the opti- 
mized (at ^ = C) VV- and PV-like algorithms is presented 
in Fig. 2 at the same set of time steps. These functions 
are plotted by curves marked as OVV and OPV, respec- 
tively. For the purpose of comparison, the functions cor- 
responding to the original VV and PV integrators are 
also drawn there (curves marked as VV and PV). Note 
that for the original integrators, the time step within each 
subset was chosen to be always twice smaller than that 
of the optimized versions (this condition is necessary to 
provide the same number of force recalculations during 
the same observation time), namely, h* = 0.005, 0.0025, 
0.00125, and 0.0005 (see subsets (a), (b), (c), and (d), 
respectively). Note also that within the original Verlet 
algorithms, the third- and higher-orders truncation un- 
certainties become too big at step sizes h* > 0.005. In 
particular, then the ratio of the total energy fluctuations 
to the fluctuations in potential energy (the standard ratio 
for estimating the precision of the calculations) appears 
to be more than a few per cent. For this reason, such 
large step sizes cannot be used in precise MD simulations 
and, thus, are not considered in the present study. 

As we see from Fig. 2, both the original (VV and PV) 
and optimized (OVV and OPV) algorithms exhibit very 
good stability properties (the excellent stability can be 
explained |3[|^ by the symplecticity and time reversibility 
of the produced solutions). No systematic deviations in 
the total energy fluctuations can be observed for all the 
integrators. Instead, in each of the cases the amplitude 
of these deviations tends to its own value which does not 
increase with further increasing the length of the simu- 
lations. However, this value appears to be significantly 
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larger for the original versions VV and PV. On the other 
hand, using the optimized OVV and OPV algorithms 
even with double sizes of the time step allows us to de- 
crease the unphysical energy fluctuations approximately 
in factor three. This is in an excellent accord with our 
theoretical prediction (14). Note also that the OPV al- 
gorithm is slightly better in energy conservation than its 
OVV version (whereas the VV integrator is better with 
respect to the PV counterpart). Furthermore, in view of 
the structure Eqs. (15) and (16), the OPV algorithm is 
more convenient when averaging macroscopic quantities. 
In particular, then the interparticle potentials can be cal- 
culated at the end of time steps simultaneously with the 
interparticle forces within the same loop, increasing the 
time efficiency of the computations. 




FIG. 2. The total energy fluctuations as functions of the 
length of the simulations performed using the optimized VV- 
(solid curve marked as OVV) and PV- (dashed curve, OPV) 
algorithms, as well as the original VV (solid curve, VV) and 
PV (dashed curve, PV) integrators. The results correspond- 
ing to different values of the time step, namely, h* — 0.01 
and 0.005, 0.005 and 0.0025, 0.0025 and 0.00125, as well as 
0.001 and 0.0005 are presented in subsets (a), (b), (c), and 
(d), respectively. 

In conclusion, we point out that new second-order 
velocity- and position- Verlet-like algorithms have been 
proposed to improve the efhciency in MD simulations of 
classical systems. The algorithms are explicit (i.e., do 
not require any iteration) , simple in implementation, and 
produce stable solutions which (like exact phase trajec- 
tories) are symplectic and time reversible. The main ad- 
vantage of the introduced algorithms with respect to the 
widely used Verlet integrators is the possibility to gener- 
ate more precise trajectories at the same overall compu- 
tational efforts. As has been demonstrated in a particular 
case of microcanonical MD simulations of a LJ fluid, the 
new algorithms allow one to reduce in several times the 
unphysical fluctuations of the total energy. 



It has been shown rigorously within a consequent the- 
oretical approach that the proposed algorithms with re- 
spect to their time efficiency should be considered as op- 
timal among all decomposition second-order integrators 
at single splitting of the Liouville operator. The opti- 
mized algorithms can be adapted to multiple scale inte- 
gration (at the presence of long-range interactions when 
the potential part of the Liouville operator is decomposed 
additionally) and extended to many-component systems 
with orientational degrees of freedom. Moreover, the pre- 
sented decomposition (8) of noncommutative operators 
is applicable for quantum Monte-Carlo simulations |T^ ] 
(because all the time coefhcients at the exponential prop- 
agators remain positive in the optimized regime) . These 
and other questions will be considered in further investi- 
gations. 
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